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Abstract 

^ In this paper, we propose a progressive Bayesian procedure, where the measurement information 

^ is continuously included into the given prior estimate (although we perform observations at 

^ discrete time steps). The key idea is to derive a system of ordinary first-order differential 

equations (ODE) by employing a new coupled density representation comprising a Gaussian 
1^ density and its Dirac Mixture approximation. The ODE is used for continuously tracking the 

^ true non-Gaussian posterior by its best matching Gaussian approximation. The performance 

^ of the new filter is evaluated in comparison with state-of-the-art filters by means of a canonical 

[ — I benchmark example, the discrete-time cubic sensor problem. 

>^ 

^ 1. Introduction 

We consider state estimation in discrete-time stochastic nonlinear dynamic systems. Thanks 
^ to their simplicity, Gaussian filters, i.e., filters representing all state densities by Gaussians, 

^ are an attractive tool for solving this type of estimation problem. However, in general, their 

^ full estimation performance is not exploited due to additional assumptions and simplifications. 

Q In this paper, we focus on Gaussian filters that operate by finding the best-matching 

^ Gaussian to the true posterior by means of an explicit shape optimization after every processing 

O step. A filter of this type is called Gaussian-assumed density filter (GADF). 

In case of the filter step, the actual measurement is considered during the approximation, 
which leads to better results than the commonly used type of Gaussian filter, the so called 
Linear Regression Kalman Filter (LRKF) [Tj. LRKFs approximate the (non-Gausian) joint 
density of state and measurement corresponding to the prior density, the noise density, and 
the given nonlinear measurement equation by a jointly Gaussian density, which is independent 
of the actual measurement and typically is only a very rough approximation of the true joint 
density. 

Examples of LRKFs |2J are the Unscented Kalman Filter (UKF) and its scaled version 
its higher-order generalization [S], a generalization to an arbitrary number of deterministic 
samples placed along the coordinate axes 0, filters performing an analytic or semi-analytic 
calculation of the required moments [7], [H] based on a decomposition into parts that can be 
calculated in closed form or via a sample approximation [0], and filters based on numerical 
integration for calculating the required nonlinear moments of the prior Gaussian density [TU] . 

Of course, GADFs are more complicated to implement in comparison to the LRKFs and 
there are various options for minimizing the shape deviation between the true posterior and 
its Gaussian approximation. One option is to employ moment matching, i.e., using the mean 
and covariance matrix of the true posterior as parameters for the desired Gaussian, as this is 



known to minimize the KuUback-Leibler distance between the two densities. Unfortunately, in 
the case of nonhnear measurement equations and the corresponding comphcated Likehhood 
function, it is in general not a simple task to calculate mean and covariance matrix of the 
true posterior, as analytic solutions are a rare exception. In most cases, numerical integration 
over the true posterior, i.e., the product of the (Gaussian) prior density and the Likelihood 
function, is in order, such as Monte Carlo integration [llj . 

In this paper, we propose a progressive Bayesian procedure for Gaussian-assumed density 
filtering, where the measurement information is continuously included into the given prior 
estimate (although we perform observations at discrete time steps). The first progressive 
filtering procedure has been developed in 2003 [|12j for state estimation with a Gaussian 
Mixture representation in the scalar case, where a homotopy continuation approach was 
proposed for tracking the true posterior with its approximation minimizing a squared-integral 
distance measure. The multi-dimensional case was treated in [13]. A generalization of this 
method to various another distance measures is proposed in [H]. Besides state estimation, the 
progressive processing idea has been applied to moment calculation [TS] and Gaussian Mixture 
reduction [161 ■ homotopy-based filtering method operating on a particle representation is 
given in |17j . 

This paper is a further development of the procedure [13] that was developed for state 
estimation with a Gaussian Mixture representation. The key idea is to derive a system 
of ordinary first-order differential equations (ODEs) by employing a new coupled density 
representation comprising a continuous and a discrete part, specifically a Gaussian density and 
its Dirac Mixture approximation [TS]. Dirac Mixture approximations for nonlinear estimation 
and filtering have been proposed for the case of scalar continuous densities in [TH], [20] An 
algorithm for sequentially increasing the number of components is given in [21] and applied to 
recursive nonlinear prediction in [22] . Mult i- dimensional Gaussian densities are treated in |18] . 
Of course, more complicated continuous densities can be handled with this approach such as 
Gaussian Mixture densities, which is outside the scope of this paper. The system of ordinary 
first-order differential equations is used for continuously tracking the true non-Gaussian 
posterior by its best matching Gaussian approximation. 

In contrast to [I3j, where a squared-integral distance measure between the true and the 
approximate posterior was minimized, here we derive ODEs for directly tracking the mean 
and covariance of the true posterior density by its Gaussian approximation. 

The new progressive estimation method allows for arbitrary noise structures, even for noise 
structures that cannot easily be treated by LRKFs such as multiplicative noise. The required 
integrals are solved by employing the discrete part, i.e., the Dirac Mixture, of the coupled 
discrete-continuous-approximation, which allows to automatically place Dirac components 
solely in the interesting regions of the state space, that is the support of the true posterior. As 
a result, the new filtering method is fast, efficient, and robust. Its performance is evaluated 
in comparison with state-of-the-art filters by means of a canonical benchmark example, the 
discrete-time cubic sensor problem. 

2. Problem Formulation 

We consider the general problem of estimating the hidden state of a discrete-time stochastic 
nonlinear dynamic system based on noisy measurements, which consists of a prediction step 
(or time update) employing a system model for propagating the estimated state from time 
step to time step and filter step (or measurement update) for including observations taken at 
a given time step into the state estimate. Here, the focus is on the filter step that is typically 
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considered harder compared to the prediction step. The insights obtained for the filter step 
can, however, be used for the prediction step as weU. 
A generative measurement equation 

y = h{x,v) (1) 

is investigated, with state x, a specific measurement y, and measurement noise v with 
corresponding noise density fv{v). The special case of additive noise 

y = h{x) + V (2) 

is also of interest, as it usually simplifies matters. 

We assume that the generative model can somehow be converted to a probabilistic model 
represented by the conditional density f{y\x). For a given specific measurement y, this 
conditional density is the so called Likelihood functioi]]^ abbreviated as 

/lU) = f{y\x) 

For the case of additive noise, the Likelihood function is given by 

fiix) = f^{y-h{x)) . 
Other noise structures result in different Likelihood functions. 

Gaussian Filters 

For a Gaussian-assumed density filter, we have a Gaussian prior density f^{x) that 
undergoes a Bayesian filter step according to the following multiplication with the Likelihood 
function 

r{x) = c-nx)-hix) , 

where the resulting posterior density is denoted by f^{x)- c is a normalization constant 
ensuring that the posterior density integrates to one. The tilde is used to underline that this 
is the true density resulting from performing a single filter step. 

Obviously, the true posterior density in general is not Gaussian anymore. In order to enable 
recursive processing without increase in computational complexity, the true posterior has to 
be approximated by a Gaussian density for the next processing step. The state-of-the-art for 
doing so will be reviewed in the next section. 

3. Density approximations 

Estimating the state of stochastic nonlinear dynamic systems poses three challenges: 1. For 
a given finite-parameter representation of the state density, the density type is usually not 
maintained across the filter or the prediction step, so that continuous re-approximations are 
necessary to allow for recursive processing. 2. By doing so, the number of parameters of the 
selected representation typically increases without bound. 3. The approximation of the true 
density by the selected representation cannot be efficiently performed or the true density 
cannot even be calculated at all. 



^The Likelihood function is not necessarily a valid density function. Although it is always non-negative, it 
does not necessarily integrate to one, i.e., it usually is not normalized or even cannot be normalized. 
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Continuous Density Representations 

For representation purposes, continuous representations such as mixtures of Gaussians are 
convenient when it comes to presenting and using the results of the state estimation process. 
They smoothly cover the considered part of the state space, easily represent complicated 
density shapes, and require relatively few parameters for satisfactory approximation results. 
However, re-approximating the true density resulting from a certain processing step by a 
continuous is usually hard. One example in the context of the Bayesian filter step is the 
re-approximation of the product of a prior Gaussian Mixture and a complicated Likelihood 
function by a posterior Gaussian Mixture. Another example concerning the prediction step is 
the propagation of a Gaussian Mixture through a nonlinear dynamic system. 

Discrete Density Representations 

The difficulties with processing continuous densities led to the development of several types 
of filters employing discrete density representations over continuous domains. Using discrete 
samples instead of a continuous representation simplifies both processing steps significantly. 
Especially the prediction step becomes almost trivial as it just consists of propagating point 
values independently through the system equation. The filter step, on the other hand, is more 
challenging. At first, it seems simple to implement as it just consists of multiplying the discrete 
samples with the Likelihood function. However, it 1. requires an explicit Likelihood function 
that is not always available and 2. in a naive implementation usually leads to a degeneration 
of the sample set as most samples obtain zero weights. 

Furthermore, it is usually difficult to assess the quality of a sample representation, which 
is related to specifying the number of samples required for a given task. Although summary 
statistics such as mean and covariance of the sample set can easily be calculated, more useful 
measures such as entropy or smoothness cannot easily be obtained. 

Although the solution to the stochastic state estimation problem consists of fully determin- 
istic densities, particle filters employing independent random samples become popular due 
to their simplicity. However, random samples of small size do not provide a homogeneous 
coverage of the underlying density. As a result, the samples usually do not adequately represent 
the underlying true density. In addition, the quality of the results depends on the specific 
realization and can only be assessed via Monte Carlo simulations and appropriate averaging. 

Alternating Density Representations 

This includes sample-based Gaussian filters, such as the Unscented Kalman Filter (UKF) 
[25] and the Gaussian filter in where a prior Gaussian density is approximated by deter- 
ministic samples that undergo the processing step. The processing result is approximated 
by a Gaussian density again and the process is repeated. This procedure cannot be easily 
generalized to more general continuous representations such as Gaussian Mixtures. This is 
due to the fact that approximating the discrete representation by a continuous representation 
is equivalent to density estimation from samples, which is a tough problem. 

New coupled discrete- continuous- approximation 

We now introduce a new class of density representations, the so called coupled discrete- 
continuous- approximations, where two types of density approximations are maintained simul- 
taneously. Instead of switching back and forth between representations convenient for the 
current task at hand, the desired representation is always available when needed. In doing 
so, we can make use of the best features of both worlds. The continuous part is employed for 
maintaining the smoothness of the representation and for taking the derivatives with respect to 
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the density parameters. The discrete part is employed for evaluating and comparing densities 
and for solving integrals in a quasi Monte-Carlo fashion. 

The principle is to have an inner core representation, a continuous one called fc, that 
controls an outer representation, a discrete one called fd- Typically, the parameters t]^ of 
the inner representation are specified and lead to appropriate parameters t]^ of the outer 
representation that is then used for, e.g., numerical integration. 

The two approximations are tied in such a way that changing the shape of one representation 
directly influences the shape of the other one and its corresponding parameters. The degrees 
of freedom of the coupled discrete-continuous-approximation are given by the number of 
parameters of the inner representation. The number of parameters of the outer approximation 
is a least equal to the number of parameters of the inner approximation, but in most cases 
will be larger. 

As a specific inner representation fdxyrjj, we will focus on Gaussian densities, where the 
parameter vector rj^ contains the elements of the mean and the covariance matrix. The outer 
representation is given by the Dirac Mixture approximation of the inner Gaussian calculated 
as described in [TS] . 

Modifying the inner parameters of the new coupled discrete-continuous-approximation 
based on employing the outer parameters for numerical integration is demonstrated in the 
next section. 



4. Progressive Filter Step 

The first step is to redefine the Likelihood function in order to achieve a continuous 
execution of the filter step. This progressive Likelihood is defined by 

where 7 is an artificial time with 7 G [0, 1]. It is desired that 
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holds. 

Several options for defining progressive Likelihood functions exist. This includes progres- 
sively modifying the given nonlinear mapping h{x) of the underlying generative model Q 
or varying the noise variance as in [13]. Here, we use the exponentiation of the Likelihood 
function as used in [T7|, [H]. The modified Likelihood function is then given by 

h{x,^) = [h{x)r . (3) 

Remark 4-1 (Additive Noise) For the special case of generative systems suffering from 
additive noise in we obtain^ 



fL{x, 7) = exp <^ - - 7 (y- h{x)) C^^ (y - h{x)) 



^For the additive noise case, exponentiation is equivalent to progressive modification of the noise covariance 
matrix as proposed in [13, since we have 

7 (y - hix)f C-i {y_ - h{x)) = {y- h{x)f f-C„) {y_ - h{x)) . 
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where multipliers have been omitted to achieve fhi^il = 0) = 1. 

The second step is the continuous execution of the filter step, for now written with a general 
prior density f^{x) and a modified posterior density /'^(x, 7) depending on the artificial time 
7 introduced above 

ffe7) = /%)-/Lfe7) 

for 7 G [0,1]. For the final tim^ 7 = 1, the modified Likelihood reaches the original 
Likelihood and as a result, the original posterior is attained by the modified posterior, i.e., 
/"^(x, 7 = 1) = /^(x). On the other extreme, for the start time 7 = 0, the modified posterior is 
desired to be identical to the prior density, i.e., we have /'^(x, 7 = 0) = f^{x). Hence, we have 

f(x,7) = ffe7 = 0)-^(x,7) 

for 7 G [0,1]. 

Plugging in the continuous part of the coupled discrete-continuous-approximation intro- 
duced in the previous Section gives 

fcix, ^ fix, 7) = Fix) ■ 7) (4) 

for 7 G [0, 1], where left-hand-side and right-hand- side now become equal only for 7 = 0, but 
only approximatelj]^ equal for 7 > 0. 

In order to find the best matching Gaussian approximation f^ix, ^7^(7)) in (jlj) we desired 
the first moments to be equal as 

m /c vS^)) dx= ^mFix) ■ 7) dx , 

with m = [l,x,x'^]'^. Taking the derivative with respect to 7 on both sides gives 

dfc{x,%{l)) f dfiix,-/) 
m -—^ ax= m j^[x) ■ ax , 



with 



We obtain 



which leads to 



(97 7r^ (97 

/Lfe7) 
pT(x,j) 1(7) 



rnp^{x,-f)dxr] h) = / m F {x) ■ f l{x, ^i) dx 

R^ - -'^ JR^ 



P(7) 

P(7) ^(7) = m ■ (5) 



•^Of course, time here denotes the artificial time 7 introduced above. 

^As long as working with the infinite-dimensional functional representation of the true density, the 
two sides in ^ are equal as /'^(a;, 7) is capable of following changes in the right-hand-side exactly. On 
the other hand, for a finite-dimensional representation f^{x,rj (7)), i.e., a density function depending on a 
finite-dimensional parameter vector rjj^pi), the left-hand-side cannot necessarily exactly follow the changes of 
the right-hand-side as the product on the right-hand-side typically is not of the same density type. 
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For the specific progression from ([3]), we obtain 



dfLix,-f) dfl{x) 



fl{x) ■ log (/^(x)) = 7) ■ log . 



Remark ^.2 (Additive Noise) A further simplification is achieved hy focusing on the case 
of additive noise. Taking the derivative of the Likelihood now gives 

hix, 7) = hix, 7) {y - h{x)Y {y - h{x)) . 

5. Scalar Gaussian- Assumed Density Filter 

This section is devoted to deriving specific formulas for scalar Gaussian densities. 

Closed-form Expressions for Matrix P(7) 

We start with a scalar Gaussian density for representing given by 



fcix,rih)) = w(7) 



1 



27r (t(7) 



exp 



1 (x — m(7))^' 

2 ^ 



with 



^,(7) = [^(7), "^(7), (^{i)y 



The required derivatives of the continuous part of the CoDiCo-approximation with 
respect to the density parameters collected in the vector p{x, 7) given by 



p(x,7) 



or 



p(a;,7) 



1 

u>(7) 
x—m{j) 



9«){7) 

8/g(^.!?,(7)) 
9m(7) 



(6) 



For this specific choice of continuous representation the matrix P(7) is given by 

/I \ 

P(7) = m(7) '"^(7) 

ym^(7) + 0-^(7) 2^(7)771(7) 2?i;(7) cr(7)y 

Solving for 17^(7) in (|5| could be performed based on the matrix P(7) directly. However, it 
is possible to give its inverse in closed form, which will be denoted by 



Q(7) = P(7)-' 



7 



with 



Q(7) 



/ 1 

m(7) 
w{-f) 

m^(7)-cr^(7) 



1 

w{-f) 



\ 




Hence, we obtain 



1 

\ 2w)(7)(7(7) w(7)s(7) 2m(7)o-(7)/ 



^(7) = Q(7)&(7) 



(7) 



with 



ujci) = [^(7), "^(7), ^^(7)]^ • 

The right-hand-side of ([T]) that we will now denote by r(7) can now be written as 

r(7) = Q(7)&(7) 

= Q(7) / m- r{,x,'-f) -log {fL{x)) dx 

q{x,l) ■ r{x,-f) -logifLix)) dx , 



R 



with 



g(x,7) 



1 

a;— 771(7) 

«!(7) 



(x-m(7))^-(a(7))^ , 
V 2«,(7)a(7) / 



Approximate Expressions for Vector 6(7) 

Calculating r(7) in ^ amounts to calculating certain nonlinear moments of the true 
posterior density, as closed-form solutions are a rare occasion, we have to be content with 
approximate integration, which is pursued further here. 

One viable option would be to replace the prior density /^(x) by its Dirac Mixture 
approximation and thus, turn the integration into summation. However, this is only efficient 
for small values of 7 as long as the progressive Likelihood function is wide and does not 
modify the prior density too much. For larger 7, the Likelihood function typically becomes 
narrower and narrower and would force many Dirac components to zero, which is similar to 
the degeneration problem in particle filtering. By doing so, only very few Dirac components 
would really contribute to the integration, especially for larger 7, which is not desired. 

Thanks to progressive processing, we have an approximate posterior fc{x,ri^{'-f)) available 
at every 7. This approximate posterior now allows to integrate only over those portions of 
the state space that contains the true posterior. For doing so, we replace the continuous part 
/^(x, 77^(7)) of the CoDiCo-approximation in (Sj) by the coupled discrete density /J, i.e., a 
scalar Dirac Mixture approximation according to 



fdi.x.%in)) = Ylwi{rih)) (5(x-Xi(r7 (7))) . 



i=l 



As a result, we obtain 



= ^Wi ■q{xi,'j) ^^^^ -logifiixi)) dx 

i=l 



/cfe ^.(7)) 



(9) 
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where the exphcit dependency of the parameters of the discrete density representation upon 
the parameters of the continuous part has been omitted. 

It is important to note that by employing Dirac components solely in regions of the 
state space where the true posterior f^{xi,'y) is represented by its Gaussian approximation 
fc{3^iVjC())i components contribute to the integration. As a result, by far fewer 

components are required to achieve a good accuracy. 

For the special case of a scalar measurement equation with additive Gaussian noise leading 
to a scalar Likelihood function 



= exp 

we obtain 



l{y-h{x)f 



2 (T, 



2 



l {y-h{ xi)f 

for 1 = 1, 



log(/L(xi))= ^ 



6. Evaluation 

For demonstrating the significant increase in performance achieved by the new filter, it is 
evaluated in comparison with the state-of-the-art filters. As a benchmark example, we consider 
the discrete-time cubic sensor problem 

y = x^ + V , 

where v is zero-mean Gaussian measurement noise. This is a canonical example, simple to 
understand, and, compared to some contrived artificial example, allows a good assessment of 
filter performance. 

For comparison purposes, the proposed new filter is compared to 

• the class of Linear Regression Kalman Filters (LRKF), specifically to 

— the Unscented Kalman Filter (UKF), 

— the Gaussian Filter (GF), 

— the Analytic LRKF, 

• and an assumed-Gaussian density filter using moment matching for fitting a Gaussian 
to the true posterior employing Monte Carlo integration for approximating the desired 
moments (MC). 

For a measurement y = 3, a prior Gaussian density with mean —1 and variance 1, and a 
cubic measurement equation corrupted by zero-mean Gaussian noise with variance 1.2, the 
results of one filter step for the different filters are shown in Figure [T] The true posterior of 
the Bayesian filter step is calculated on a fine grid with 30.000 grid points. The ground truth 
for comparison purposes is the Gaussian density with mean and variance identical to the true 
posterior. A standard Unscented Kalman Filter (UKF) is employed. The Gaussian Filter 
(GF) is employed based on an optimal Dirac Mixture approximation of the prior Gaussian 
density according to [6J, where the filter step is performed by assuming that measurement and 
state are jointly Gaussian. The Analytic LRKF computes the moments in closed form and 
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Figure 1: Comparison of the results for one filter step for the ground truth, i.e., the best Gaussian approximation 
computed numerically, compared with the true posterior, the Gaussian Filter (GF), the Unscented Kalman 
Filter (UKF), the Analytic LRKF, the Gaussian Particle Filter (GPF), and the new progressive Gaxissian filter. 



additionally assumes that measurement and state are jointly Gaussian. The Gaussian Particle 
Filter (GPF) computes posterior mean and variance by means of 10^ random samples drawn 
from the prior Gaussian density that are used for evaluating the product of Likelihood and 
prior density. Subsequently, the resulting weighted samples are used for moment matching. 
The new progressive Gaussian filter is applied with a discrete density approximation using 30 
Dirac components. 

Figure [T] clearly shows that the LRKFs have difficulties approximating the true posterior. 
In contrast, the GPF for 10^ samples and the new progressive Gaussian filter produce results 
almost identical to the ground truth. 

To illustrate the savings in the number of samples used by the new progressive Gaussian 
filter, its estimation quality is compared to the Gaussian particle filter in Figure [2] for the 
same experimental parameters as above. It is obvious that very good estimation results are 
already obtained for very few (deterministic) samples. In contrast, the results of the Gaussian 
particle filter are non-deterministic and the average convergence is much slower. 

The results of recursive filtering over 50 times steps are shown in Figure [3] The recursion 
is started with a prior Gaussian density with mean —1 and variance 30. Noise mean is 0, the 
noise variance is 1.2. From time step 1 to 19, the true state is 1. At time step 20, the true 
state is changed to and Gaussian noise with variance of 9 is added to all estimates. 

The top plot in Figure |3] shows the estimated means of the true posterior, its best Gaussian 
approximation^ the new progressive Gaussian filter, and the Analytic LRKF. The other 
variants of LRKF have been omitted as the Analytic LRKF provides their lower bound in 
terms of estimation quality. 



^The best Gaussian approximation used as the ground truth at every time step is also recursively calculated 
on the grid based on the previous best Gaussian approximation. The accumulating error between the true 
posterior and the ground truth is unavoidable due to the Gaussian assumption. 
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Figure 2: Comparison of the estimation quality as a function of the number of samples used between the 
Gaussian particle filter (GPF) and the proposed progressive Gaussian filter. As the GPF provides non- 
deterministic results, 100 Monte-Carlo runs have been performed to produce the RMSE and the min-max error 
bars. 

The bottom plot in Figure |3] shows the absolute mean error of the best Gaussian approxi- 
mation, the new progressive filter, and the Analytic LRKF with respect to the mean of the 
true posterior. 

The two plots in Figure [3] show that the proposed progressive Gaussian filter provides 
results very close to the best Gaussian approximation, which itself is rather close to the mean 
of the true posterior. The analytic LRKFs represented by their highest-quality variant, the 
Analytic LRKF, show a much larger deviation. 

7. Conclusions 

A progressive Gaussian filtering method has been introduced that is based on a coupled 
discrete-continuous-representation with a Gaussian density as the continuous part and its Dirac 
Mixture approximation as the discrete part. For the moment calculation. The Dirac Mixture 
part is employed to evaluate the true posterior density at discrete points only. In contrast to 
Monte-Carlo methods available for that purpose, the Dirac components all contribute to the 
integration as progressive processing is exploited to solely place components in the important 
regions of the state space. As a result, the new filter requires significantly less samples, is 
fast, efficient, and robust, and achieves a high estimation quality compared to state-of-the-art 
Gaussian filters. It can easily be used as a replacement for standard Gaussian filters such as 
the UKF. 
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